## This code imports and process data.

# Import libraries ----

library(raster) # Process raster data
library(rgdal) # Import raster and vector data
library(rgeos) # Topology operations
library(sp) # Classes and methods for spatial data
library(conleyreg) # Conley standard errors
library(geosphere) # Measure geographic distance
library(rmapshaper) # Simplify shapefiles
library(robustHD) # Standardize variables
library(readxl) # Read excel files
library(stargazer) # Latex table
library(varhandle) # Call variables for loops
library(xtable) # Create result tables
library(ggmap) # Import Google Maps
library(ggplot2) # Create plot

# Import data ----

# ** Secondary data ----

# **** Import 

konbaung <- ms_simplify(readOGR("konbaung","konbaung_1817")) # Konbaung Dynasty porders

myanmar <- ms_simplify(readOGR("mmr_polbnda2_adm1_250k_mimu","mmr_polbnda2_adm1_250k_mimu")) # Bago village tracts

altitude <- raster("fao_geo/fao_median_altitude/lr_ter_altmed05.tif") # Median altitude (in meters)

slope_idx <- raster("fao_geo/fao_slope_index/lr_ter_slpidx05.tif") # Slope index as a measure of terrain ruggedness (0-1 scale)

soil_quality <- raster("fao_geo/fao_soil_quality/lr_soi_sq7b_mze.tif") # Soil quality (ordinal,larger numbers mean lower quality)

rivers <- ms_simplify(readOGR("myanmar_river_network_250k","myanmar_river_network_250k")) # Major rivers

# **** Process

wgs84 <- crs(myanmar) # Extract WGS 1984 projection

region <- myanmar[myanmar$ST %in% c("Bago (East)","Bago (West)"),]

region <- gUnionCascaded(region)

myanmar <- gUnionCascaded(myanmar)

altitude <- intersect(altitude,region) # Altitude raster for Bago region

slope_idx <- intersect(slope_idx,region) # Slope raster for Bago region

soil_quality <- intersect(soil_quality,region) # Soil quality raster for Bago region

# ** Original 1784 sittan data ----

# **** Import 

sittan_towns <- read.csv("sittan_towns_1783.csv") # Town data from 1784 sittan

sittan_towns <- sittan_towns[complete.cases(sittan_towns$hereditary),]

sittan_pops <- read.csv("sittan_pop_1783.csv") # Population data from 1784 sittan

sittan_pops <- sittan_pops[complete.cases(sittan_pops$longitude),] # Include only towns with available geo-location

# **** Process 

sittan_towns_shp <- SpatialPointsDataFrame(sittan_towns[,c("longitude","latitude")],sittan_towns,proj4string=wgs84) #Transform to spatial vector data frame

sittan_pops_shp <- SpatialPointsDataFrame(sittan_pops[,c("longitude","latitude")],sittan_pops,proj4string=wgs84) #Transform to spatial vector data frame

sittan_towns_shp <- intersect(sittan_towns_shp,region)

sittan_pops_shp <- intersect(sittan_pops_shp,region)

# ** Original 1912 gazetteer data ----

# **** Import

colonial_towns <- read_excel("colonial_towns.xlsx") # Towns and villages from 1912 British Burma Gazetteer

colonial_towns <- colonial_towns[is.na(colonial_towns$longitude)==FALSE,] # Remove NA values

colonial_towns <- colonial_towns[is.na(colonial_towns$houses)==FALSE,] # Remove NA values

colonial_towns$coercive <- (colonial_towns$m_police + colonial_towns$c_police + colonial_towns$outpost)

colonial_towns$welfare <- (colonial_towns$school + colonial_towns$hospital)

# **** Process

colonial_towns_shp <- SpatialPointsDataFrame(colonial_towns[,c("longitude","latitude")],colonial_towns,proj4string=wgs84) #Transform to spatial vector data frame

colonial_towns_shp <- intersect(colonial_towns_shp,region) # Select towns in Bago region

# ** Original conflict data ----

# **** Import

conflict <- read.csv("conflict_only.csv") # Armed conflict events from 1948 to 1998

# **** Process

conflict <- conflict[complete.cases(conflict$long),] # Select only events with available geo-location

conflict_shp <- SpatialPointsDataFrame(conflict[,c("long","lat")],conflict,proj4string=wgs84) #Transform to spatial vector data frame

conflict_shp <- intersect(conflict_shp,region) # Conflict for Bago region

# ** Location of Yangon ----

# **** Import

Yangon <- c(96.158840,16.774383) # Location of Yangon (Sule Pagoda)

# **** Process

Yangon <- SpatialPoints(rbind(Yangon)) # Convert to SpatialPoints class

# Calculate data ----

# ** Calculate geographic covariate values for towns and villages ----

colonial_towns_shp$altitude <- extract(altitude,colonial_towns_shp,fun=mean) # Altitude

colonial_towns_shp$slope_idx <- extract(slope_idx,colonial_towns_shp,fun=mean)/10000 # Slope index

colonial_towns_shp$soil_quality <- extract(soil_quality,colonial_towns_shp,fun=mean) # Soil quality

colonial_towns_shp$soil_quality <- ifelse(is.na(colonial_towns_shp$soil_quality)==TRUE,0,colonial_towns_shp$soil_quality)

river_dist <- dist2Line(colonial_towns_shp,rivers) # Distance from a major river

colonial_towns_shp$river_dist <- river_dist[,1]/1000 #Attach distance from river

yangon_dist <- pointDistance(colonial_towns_shp,Yangon,lonlat = TRUE) # Distance from Yangon

colonial_towns_shp$yangon_dist <- yangon_dist/1000 #Attach distance from Yangon

neighbor_dist <- pointDistance(colonial_towns_shp,sittan_towns_shp,lonlat = TRUE) # Distance from the towns without appointed headmen

colonial_towns_shp$neighbor_id <- apply(neighbor_dist,1,which.min) # Attach the minimum distance from the towns without appointed headmen

colonial_towns_shp$neighbor_dist <- apply(neighbor_dist,1,FUN = function(x) min(x)/1000) # Attach the minimum distance from the towns without appointed headmen

for(i in 1:nrow(colonial_towns_shp)){
  colonial_towns_shp$replaced[i] <- sittan_towns_shp$hereditary[colonial_towns_shp$neighbor_id[i]]
  colonial_towns_shp$region[i] <- sittan_towns_shp$region[colonial_towns_shp$neighbor_id[i]]
}

colonial_towns_shp$replaced <- 1 - colonial_towns_shp$replaced

# **** Calculate post-colonial variables ----

# Count the number of armed conflicts within 10 kilometers

for(i in 1:nrow(colonial_towns_shp)){
  distance_matrix <- t(distm(colonial_towns_shp[i,],conflict_shp,fun=distGeo))
  distance_matrix <- data.frame(distance_matrix)
  colonial_towns_shp$conflict_count[i] <- length(distance_matrix[distance_matrix<=1000])
}

# Standardize variables

colonial_towns_shp$houses_sd <- standardize(log(colonial_towns_shp$houses), centerFun = mean, scaleFun = sd)

colonial_towns_shp$altitude_sd <- standardize(colonial_towns_shp$altitude, centerFun = mean, scaleFun = sd)

colonial_towns_shp$slope_idx_sd <- standardize(colonial_towns_shp$slope_idx, centerFun = mean, scaleFun = sd)

colonial_towns_shp$soil_quality_sd <- standardize(colonial_towns_shp$soil_quality, centerFun = mean)

colonial_towns_shp$yangon_dist_sd <- standardize(log(colonial_towns_shp$yangon_dist), centerFun = mean, scaleFun = sd)

colonial_towns_shp$river_dist_sd <- standardize(log(colonial_towns_shp$river_dist), centerFun = mean, scaleFun = sd)

colonial_towns_shp$coercive_sd <- standardize(colonial_towns_shp$coercive, centerFun = mean, scaleFun = sd)

colonial_towns_shp$m_police_sd <- standardize(colonial_towns_shp$m_police, centerFun = mean, scaleFun = sd)

colonial_towns_shp$total_force <- colonial_towns_shp$c_force + colonial_towns_shp$m_force

colonial_towns_shp$total_force_sd <- standardize(colonial_towns_shp$total_force, centerFun = mean, scaleFun = sd)

colonial_towns_shp$m_force_sd <- standardize(colonial_towns_shp$m_force, centerFun = mean, scaleFun = sd)

# Transform by inverse hyperbolic sine function

colonial_towns_shp$coercive_ih <- log(colonial_towns_shp$coercive + sqrt(colonial_towns_shp$coercive ^ 2 + 1))

colonial_towns_shp$m_police_ih <- log(colonial_towns_shp$m_police + sqrt(colonial_towns_shp$m_police ^ 2 + 1))

colonial_towns_shp$total_force_ih <- log(colonial_towns_shp$total_force + sqrt(colonial_towns_shp$total_force ^ 2 + 1))

colonial_towns_shp$m_force_ih <- log(colonial_towns_shp$m_force + sqrt(colonial_towns_shp$m_force ^ 2 + 1))

colonial_towns_shp$longitude_2 <- colonial_towns_shp$longitude

colonial_towns_shp$latitude_2 <- colonial_towns_shp$latitude

